*******************
* Load in dataset *
*******************

* Set filepath of the .dta file, e.g.
use "C:\~\FPA Replication\FPA_ReplicationData.dta", clear

keep approval econ_improve Year Month TotalAircrafts monthnum modate Gr_milex Tu_milex year trk_gdppc event event_kurds ln_tradevol

****************************
* Further Data Preparation *
****************************
tsset modate

* Interpolate small gaps in approval ratings and econ perceptions (three or less months)
gen approval_ip=approval
replace approval_ip=(l.approval+f.approval)/2 if approval==. & l.approval!=. & f.approval!=.

replace approval_ip=(l2.approval+1.5*f.approval)/2.5 if approval==. & l.approval==. & l2.approval!=. & f.approval!=.
replace approval_ip=(1.5*l.approval+f2.approval)/2.5 if approval==. & f.approval==. & l.approval!=. & f2.approval!=.

replace approval_ip=(l2.approval+f2.approval)/2 if approval==. & l.approval==. & f.approval==. & l2.approval!=. & f2.approval!=.
replace approval_ip=(l3.approval+2*f.approval)/3 if approval==. & l.approval==. & l2.approval==. & l3.approval!=. & f.approval!=.
replace approval_ip=(2*l.approval+f3.approval)/3 if approval==. & f.approval==. & f2.approval==. & l.approval!=. & f3.approval!=.

gen econ_improve_ip=econ_improve
replace econ_improve_ip=(l.econ_improve+f.econ_improve)/2 if econ_improve==. & l.econ_improve!=. & f.econ_improve!=.

replace econ_improve_ip=(l2.econ_improve+1.5*f.econ_improve)/2.5 if econ_improve==. & l.econ_improve==. & l2.econ_improve!=. & f.econ_improve!=.
replace econ_improve_ip=(1.5*l.econ_improve+f2.econ_improve)/2.5 if econ_improve==. & f.econ_improve==. & l.econ_improve!=. & f2.econ_improve!=.

replace econ_improve_ip=(l2.econ_improve+f2.econ_improve)/2 if econ_improve==. & l.econ_improve==. & f.econ_improve==. & l2.econ_improve!=. & f2.econ_improve!=.
replace econ_improve_ip=(l3.econ_improve+2*f.econ_improve)/3 if econ_improve==. & l.econ_improve==. & l2.econ_improve==. & l3.econ_improve!=. & f.econ_improve!=.
replace econ_improve_ip=(2*l.econ_improve+f3.econ_improve)/3 if econ_improve==. & f.econ_improve==. & f2.econ_improve==. & l.econ_improve!=. & f3.econ_improve!=.

* log Military Expenditures, GDP pc, and Violation counts
gen ln_tumilex=ln(Tu_milex)
gen ln_grmilex=ln(Gr_milex)
gen ln_tugdppc=ln(trk_gdppc)

gen ln_airvio=ln(1+TotalAircrafts)
gen airvio=TotalAircrafts

* Generate a dummy for Onset of Natural Gas Dispute
gen gas_disp=0
replace gas_disp=1 if modate>tm(2018m10)

* Dummy indicators for domestic political events.
* For turkish national elections, als includes constituional referendum
gen tu_elec=0
replace tu_elec=1 if modate==tm(2015m6) | modate==tm(2015m11) | modate==tm(2018m6) | ///
modate==tm(2014m8) | modate==tm(2017m4)

gen f_tu_elec=f.tu_elec
replace f_tu_elec=0 if f_tu_elec==.

gen tu_coup=0
replace tu_coup=1 if modate==tm(2016m7)

gen tu_gezi=0
replace tu_gezi=1 if modate>=tm(2013m5) & modate<=tm(2013m8)

gen tu_turmoil=0
replace tu_turmoil=1 if modate>=tm(2013m5) & modate<=tm(2013m8)
replace tu_turmoil=1 if modate==tm(2016m7)

gen admin_no=.
replace admin_no=1 if modate>=tm(2003m3) & modate<=tm(2007m8)
replace admin_no=2 if modate>=tm(2007m9) & modate<=tm(2011m6)
replace admin_no=3 if modate>=tm(2011m7) & modate<=tm(2014m8)
replace admin_no=4 if modate>=tm(2014m9) & modate<=tm(2018m6)
replace admin_no=5 if modate>tm(2018m6) & modate!=.
tab admin_no, gen(admin_fe)

* For greek national elections, als includes EU referendum
gen gr_elec=0
replace gr_elec=1 if modate==tm(2010m2) | modate==tm(2015m2) | modate==tm(2020m1) | ///
modate==tm(2012m5) | modate==tm(2012m6) | modate==tm(2009m10) | modate==tm(2015m1) | ///
modate==tm(2015m9)| modate==tm(2019m7)

gen f_gr_elec=f.gr_elec
replace f_gr_elec=0 if f_gr_elec==.

gen gr_euref=0
replace gr_euref=1 if modate==tm(2015m7)
gen f_gr_euref=f.gr_euref
replace f_gr_euref=0 if f_gr_euref==.

label variable airvio "Airspace Violations"
label variable approval_ip "Erdogan Approval"
label variable econ_improve_ip "Econ perception: positive"

* Seasonality variables
gen may=0
replace may=1 if monthnum==5

gen covid=0
replace covid=1 if modate>tm(2020m2)


*********************************************
* Time-series Diagnostics (see Appendix, A1) *
*********************************************
* Table A1
dfuller ln_airvio if modate>tm(2013m8)
pperron ln_airvio if modate>tm(2013m8)

* Figure A1
bysort monthnum: egen volmean=mean(airvio)
twoway (line airvio monthnum if Year==2009) (line airvio monthnum if Year==2010) ///
(line airvio monthnum if Year==2011) (line airvio monthnum if Year==2012) ///
(line airvio monthnum if Year==2013) (line airvio monthnum if Year==2014) ///
(line airvio monthnum if Year==2015) (line airvio monthnum if Year==2016) ///
(line airvio monthnum if Year==2017) (line airvio monthnum if Year==2018) ///
(line airvio monthnum if Year==2019) (line airvio monthnum if Year==2020) ///
(line volmean monthnum, lcolor(red)), legend(off) xscale(range(12)) scheme(plotplain) ///
xtitle("Month") ytitle("Airspace Violations") xlabel(1(1)12, nogrid) ylab(, nogrid)

* Table A2
tsset modate

quietly arima ln_airvio l.approval_ip f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1 2) ma(1)
estat ic
* AR2, MA1, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1 2) ma(1)
estat ic
* AR3, MA1, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1 2 3) ma(1)
estat ic
* AR3, MA2, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1 2 3) ma(1 2)
estat ic
* AR2, MA2, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1 2) ma(1 2)
estat ic
* AR2, MA3, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1 2) ma(1 2 3)
estat ic
* AR1, MA1, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1) ma(1)
estat ic
* AR1, MA2, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1) ma(1 2)
estat ic
* AR1, MA3, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1) ma(1 2 3)
estat ic
* AR3, MA3, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1 2 3) ma(1 2 3)
estat ic
* AR1, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1)
estat ic
* AR2, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1 2)
estat ic
* AR3, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ar(1 2 3)
estat ic
* MA1, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ma(1)
estat ic
* MA2, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ma(1 2)
estat ic
* MA3, May Dummy 
quietly arima ln_airvio l.approval_ip may f_tu_elec gas_disp tu_coup admin_fe4 admin_fe5 if modate>tm(2013m8), ma(1 2 3)
estat ic


*****************
* ARMA Analysis *
*****************

* Figure 1
tsset modate

twoway (line airvio modate) (line approval_ip modate, yaxis(2)) ///
if modate>tm(2013m8), scheme(plotplain) legend(position(6)) xtitle("Month") ///
ytitle("Airspace Violations", axis(1)) ytitle("Erdogan Approval", axis(2)) ///
xlabel(#7, nogrid) ylab(, nogrid)

* Table 1
sum ln_airvio l.approval_ip f_tu_elec tu_coup gas_disp l.econ_improve_ip  l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8)

* Table 2 (+ Appendix Table A3)
arima ln_airvio l.approval_ip admin_fe4 admin_fe5 may if modate>tm(2013m8), ar(1)
estat aroots, nograph
arima ln_airvio l.approval_ip admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp  if modate>tm(2013m8), ar(1)
estat aroots, nograph
arima ln_airvio l.approval_ip l.econ_improve_ip admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp if modate>tm(2013m8), ar(1)
estat aroots, nograph
arima ln_airvio l.approval_ip admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)
estat aroots, nograph


* Figure 2
arima ln_airvio l.approval_ip admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)
margins, dydx(l.approval_ip) level(90)
margins, at(l.approval_ip=(37.5(2.5)67.5)) level(90) 
marginsplot, scheme(plotplain) ciopts(recast(rline) lpattern(dash)) ytitle(ln Airspace Violations) xtitle(Approval Rate) ///
legend(off) title(" ") text(5 38.2 "{it:141}") text(4.27 68.2 "{it:68}") text(4.68 50.7 "{it:104}") 

*********************************
* VAR Diagnstotics and Analysis *
*********************************

* Table A4
dfuller approval_ip
pperron approval_ip

* Table A5
varsoc ln_airvio approval_ip if modate>tm(2013m8), exog(admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc)

* Table A6 and A7
matrix A1 = (1,0\ .,1)
matrix B1 = (.,0 \ 0,.)

svar ln_airvio approval_ip if modate>tm(2013m8), exog(admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc) ///
lags(1) aeq(A1) beq(B1) var

vargranger
varstable

* Figure 3
cd "C:\Users\mariu\Documents"
irf create order1, set(var2.irf) replace step(12)

irf graph sirf, xlabel(0(2)12) irf(order1) impulse(approval_ip) response(ln_airvio) ///
yline(0,lcolor(black)) scheme(plotplain) level(90) byopts(note("") legend(off)) ///
subtitle(" ", fcolor(none) lstyle(none)) xtitle("Months after Approval Shock")  ///
name(one, replace) ytitle("Airspace Violations Response")
irf graph sirf, xlabel(0(2)12) irf(order1) impulse(ln_airvio) response(approval_ip) ///
yline(0,lcolor(black)) scheme(plotplain) level(90) byopts(note("") legend(off)) ///
subtitle(" ", fcolor(none) lstyle(none)) xtitle("Months after Airspace Violations Shock")  ///
name(two, replace) ytitle("Approval Response")
graph combine one two, ycommon scheme(plotplain)

*********************
* Robustness Checks *
*********************

* Table A8 - No interpolation in Approval
arima ln_airvio l.approval admin_fe4 admin_fe5 may if modate>tm(2013m8), ar(1)
arima ln_airvio l.approval admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp if modate>tm(2013m8), ar(1)
arima ln_airvio l.approval l.econ_improve_ip admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp if modate>tm(2013m8), ar(1)
arima ln_airvio l.approval admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)

* Table A9 - * Coup control and structural break
list modate if tu_coup==1
gen postcoupat=0
replace postcoupat=1 if modate>=tm(2016m7)
gen appXcoup=l.approval_ip*postcoupat
arima ln_airvio l.approval_ip postcoupat admin_fe4 admin_fe5 may f_tu_elec gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)
arima ln_airvio l.approval_ip postcoupat appXcoup admin_fe4 admin_fe5 may f_tu_elec gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)

* Tables A10+A11 - Further political events
gen deadlyearthquake=0
replace deadlyearthquake=1 if modate==tm(2017m7) | modate==tm(2017m6) | modate==tm(2015m11) ///
| modate==tm(2014m5) | modate==tm(2020m6) | modate==tm(2020m2) | modate==tm(2020m1) | modate==tm(2019m9)
gen visit = 0
replace visit=1 if modate==tm(2017m12) | modate==tm(2019m2)
gen davut_exit=0
replace davut_exit=1 if modate>=tm(2016m5)
gen natosummit=0
replace natosummit=1 if modate==tm(2010m11) | modate==tm(2012m5) | modate==tm(2014m9) | modate==tm(2016m7) ///
| modate==tm(2017m5) | modate==tm(2018m7) | modate==tm(2019m12)
gen syria_op=0
replace syria_op=1 if modate==tm(2019m11) | modate==tm(2019m10)
replace syria_op=1 if modate>=tm(2018m1) & modate<=tm(2018m3)
replace syria_op=1 if modate>=tm(2016m8) & modate<=tm(2017m3)

arima ln_airvio l.approval_ip f_gr_elec admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)
arima ln_airvio l.approval_ip deadlyearthquake admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)
arima ln_airvio l.approval_ip visit admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)
arima ln_airvio l.approval_ip davut_exit admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)

arima ln_airvio l.approval_ip natosummit admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)
arima ln_airvio l.approval_ip syria_op admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)
arima ln_airvio l.approval_ip l.event admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)
arima ln_airvio l.approval_ip l.event_kurds admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)
arima ln_airvio l.approval_ip l.ln_trade admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc if modate>tm(2013m8), ar(1)

* Table A12 - AR Poisson
* If arpois package is not installed: net install sts13.pkg

gen lag_approval=l.approval_ip
gen lag_econperc=l.econ_improve_ip
gen lag12_tumil=l.l12.ln_tumilex
gen lag12_gremil=l.l12.ln_grmilex
gen lag12_gdp=l.l12.ln_tugdppc

arpois airvio lag_approval admin_fe4 admin_fe5 may if modate>tm(2013m8), ar(1) delete
arpois airvio lag_approval admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp if modate>tm(2013m8), ar(1) delete
arpois airvio lag_approval admin_fe5 may f_tu_elec tu_coup gas_disp lag_econperc if modate>tm(2013m8), ar(1) delete
arpois airvio lag_approval admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp lag12_tumil lag12_gremil lag12_gdp if modate>tm(2013m8), ar(1) delete

* Table A13 - Lag Order Four + reversed contemporaneous effects
svar ln_airvio approval_ip if modate>tm(2013m8), exog(admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc) ///
lags(1 2 3 4) aeq(A1) beq(B1) var
vargranger
svar  approval_ip ln_airvio if modate>tm(2013m8), exog(admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc) ///
lags(1) aeq(A1) beq(B1) var
vargranger

* Figure A2
svar ln_airvio approval_ip if modate>tm(2013m8), exog(admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc) ///
lags(1 2 3 4) aeq(A1) beq(B1) var
irf create order1, set(var2.irf) replace step(12)
irf graph sirf, xlabel(0(2)12) irf(order1) impulse(approval_ip) response(ln_airvio) ///
yline(0,lcolor(black)) scheme(plotplain) level(90) byopts(note("") legend(off)) ///
subtitle(" ", fcolor(none) lstyle(none)) xtitle("Months after Approval Shock")  ///
name(one, replace) ytitle("Airspace Violations Response")
irf graph sirf, xlabel(0(2)12) irf(order1) impulse(ln_airvio) response(approval_ip) ///
yline(0,lcolor(black)) scheme(plotplain) level(90) byopts(note("") legend(off)) ///
subtitle(" ", fcolor(none) lstyle(none)) xtitle("Months after Airspace Violations Shock")  ///
name(two, replace) ytitle("Approval Response")
svar  approval_ip ln_airvio if modate>tm(2013m8), exog(admin_fe4 admin_fe5 may f_tu_elec tu_coup gas_disp l12.ln_tumilex l12.ln_grmilex l12.ln_tugdppc) ///
lags(1) aeq(A1) beq(B1) var
irf create order1, set(var2.irf) replace step(12)
irf graph sirf, xlabel(0(2)12) irf(order1) impulse(approval_ip) response(ln_airvio) ///
yline(0,lcolor(black)) scheme(plotplain) level(90) byopts(note("") legend(off)) ///
subtitle(" ", fcolor(none) lstyle(none)) xtitle("Months after Approval Shock")  ///
name(three, replace) ytitle("Airspace Violations Response")
irf graph sirf, xlabel(0(2)12) irf(order1) impulse(ln_airvio) response(approval_ip) ///
yline(0,lcolor(black)) scheme(plotplain) level(90) byopts(note("") legend(off)) ///
subtitle(" ", fcolor(none) lstyle(none)) xtitle("Months after Airspace Violations Shock")  ///
name(four, replace) ytitle("Approval Response")
graph combine one two three four, ycommon scheme(plotplain)





